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Abstract 

Commonly used one-dimensional geoid models predict that the 
isostatlc geoid anomaly over old ocean basins for the boundary-layer 
therimil model of the lithosphere is a factor of two greater than that 
for the plate model. Calculations presented here, using the spherical 
analogues of the plate and boundary-layer thermal models, show that for 
the actual global distribution of plate ages, one-dimensional models are 
not accurate and a spherical, fully three-dimensional treatment is 
necessary. Tlie maximum difference in geoid heights predicted for the 
two models is only about two meters. The thermal structure of old 
lithosphere is unlikely to be resolvable using global geoid anomalies. 
Stripping the effects of plate aging and a hypothetical uniform, 35 kra, 
isostat ical ly-compensated continental crust from the observed geoid 
emphasizes that the largest -ampl i tude geoid anonuily is the geoid low of 
almost 120 m over West Antarctica, a factor of two greater than the low 


of 60 ra over Ceylon. 
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1* Introduction 

One of the ntore striking successes of the kinematic model of 
plate tectonics is that the variation in oceanic depth away from ridges 
can be understood in terms of simple models of cooling of the 
lithosphere (I'S)* Oceanic heat flow can be explained by these models 
if hydrothermal activity at the ridge is taken into account (6); 
recently it has been demonstrated that the short wavelength change in 
geoid elevation over ridges predicted by these simple models is also 
readily observed (7-9). 

A number of models have been proposed for the thermal 
evolution of the lithosphere* In the boundary layer model, the 
uppermost mantle is supposed to cool by conduction, with the 
lithospheric thickness increasing as the square root of age (2,5)* In 
the plate model, the temper-'iture is assumed to be held fixed by some 
mechanism at a depth on the order of 100 km (1,3,4). The models are 
quite similar near the ridge crest before cooling has penetrated very 
deep, but over the ocean brsins, the models differ substantially* 

Other thermal models have been proposed which Include the 
effects of latent heat from freezing partial melts (10), pressure 
gradients from return flow (11,12), heat flux into the base of the 
lithosphere (13), and heat sources at the base of or below the 
lithosphere (14,15)* Any or all of these effects may be Important in 
the earth, but in this paper, discussion will be focused primarily on 
the differences between the plate and boundary layer models* 

Three types of observations have been proposed to provide 
possible constraints on the thermal structure of the lithosphere* Heat 
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flow Is proportional to the surface temperature gradient* The surface 
temperature gradient* for the plate and boundary layer models differ 
only slightly for lithosphere less than 200 Myt old, however, so heat 
flow Is not at present useful for distinguishing between these models 
(16). 

Seafloor elevation Is related to the average heat content of 
the underlying column of rock* As the uppermost mantle cools while 
moving away from the ridge, It contracts, leading to an Increase In 
ocean depth. The total Increase In ocean depth with time (assuming 
constant thermal expansivity) Is proportional to the change In the 
average heat content of that part of the mantle which Is dynamically 
coupled to the surface, not. In general, the lithosphere alone (IS). 
This change In heat content Is Just the difference between the 
time-integrated heat flux out the system and any Internal heating or 
flux of heat Into the system* The standard boundary layer and plate 
models (assuming dynamic decoupling at the base of the plate) predict 
different seafloor elevations for old lithosphere; although the heat 
flow out of each Is aearly Identical, the plate model supposes that 
extra heat Is added to the system at the base of the plate* This 
specifically-located heat source would be Indistinguishable from any 
more general dlstrlbialon of Internal heat sources, however, since 
topography gives only a measure of the average heat content* 

A third potentially useful observation Is geold height* In 
the one-dlmenslonal (long wavelength) limit, the geold height Is 
proportional to the first moment of the density distribution (7-9). 
Because It Is sensitive to the distribution of density, geold height In 
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principle might be expected to provide a more sensitive constraint on 
thermal models than heat flow or seafloor topography, which depend upon 
surface gradient and average heat content respectively. 

In the one-dimensional limit, the geoid heights for the two 
models discussed here differ by a factor of two over old ocean basins 
(17). The actual distribution of lithospheric ages is distinctly 
three-dimensional, however, so the validity of the one-dimensional 
calculation is not a priori p^iar ' nteed. Model geolds computed using an 
estimate of the actual distribution of aeafloor ages are presented 
below. Unfortunately, the full three-dimensional, spherical calculation 
shows that the difference between the geoids computed for the two 
llthosperic models is very small and unlikely to be observable. 

2. Calculation of Isostatic Geoid Anomalies 

Dahlen (18) has recently discussed the application of the 
concept of isostacy to a sphere, noting that the theoretical geoid 
anomaly is very sensitive to the exact definition used for isostacy. 
Isostatlc geoid anomalies are by nature first order quantities; if they 
are to be computed to first order, isostacy must be defined to second 
order. Dahlen showed that geoid anomalies could differ by more than a 
factor of two, depending on the definition chosen, e.g. , constant 
pressure at the level of compensation versus constant mass above the 
level of compensation. 

In this paper, the isostatic condition is assumed to 
correspond to a local compensation such as would be achieved by a set of 
vertical (or conical) piles of different lengths and densities decoupled 
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from each other floating in an Inviscld fluid* The tfelght of any pile 
is taken to be equal to the weight of fluid displaced by the pile* 
Self-compression in the earth causes the gravitational acceleration g to 
be essentially constant with depth* so the constant weight criterion is 
equivalent to a constant mass criterion* For a halfspace* Isostacy is 
achieved if 


j Ap dz ■ 0 ( 1 ) 

where z is depth, and Ap is the density contrast between the pile and 
the fluid* The Integral is taken over the depth of the pile, l*e*, the 
region in which density contrasts are nonzero, which, in this paper, 
corresponds to the thickness of the lithosphere, h* 

For a sphere, the equivalent definition is that 

J dz - 0 (2) 

where a is the radius of the earth* Since h « a and Ap is a first 
order quantity, (1) and (2) are equivalent to first order, but they 
result in different gravitational potential anomalies (18)* 

In calculating the geoid anomaly from an arbitrary density 
distribution, it is convenient to separate the horizontal and depth 
dependences of density* The horizontal dependence is expressed in terms 
of harmonic functions. For a flat earth, the contribution of density 
contrasts Ap|^(z), of wavelength k, extending over depth h, to the geoid 
height, Nj^(x), can be written: 
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Here x is the horizontal position vector, G the gravitational constant, 
g the gravitational acceleration, and a phase angle* 

Taking (1) as the definition of isostacy and assuming kh « 1 

leads to the simplification 
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This approximation is in error by less than 10% for kd < .2, or for 
wavelengths greater than 3800. km if h corresponds to the plate thickness 
of 120 km. It overestimates the geold height contribution from shorter 
wavelengths. 

For a sphere, horizontal density contrasts can be expressed in 
terms of surface spherical harmonics of degree i. and order m. 
The equivalent of (3) is then 
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where a is the radius of the earth* For h/a « 1, using (2) leads to 




In this case, the approximation is in error by less than lOZ for i < 12 
if h ~ 120 km* For i > 20, (6) overestimates the contribution of a 
density contrast at a depth of 120 km by 20Z* 

An additional complication arises for isostatic compensation 
on a spherical earth. The interior of the earth is assumed to be 
Inviscid on the timescales applicable to isostatic equilibrium. The 
gravitational potential from an isostatically compensated mass will 
cause a net deformation of the (fluid) earth* This deformation, in 
turn, causes an additional gravitational potential, which leads to yet 
more deformation. 

This is a classical problem in geophysics which has often been 
solved for the degree two deformation associated with earth rotation and 
tides using Love numbers (19). The total geoid anomaly Ntm is given by 

"i. ■ »» "i. 

where is the Love number for degree i* 

Exact computation of the Love numbers for a fluid sphere which 
has the density and compressibility distribution of the earth is a 
formidable task. To my knowledge, it has been done only for degree two. 
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for which H 2 ■ 1*96 (19). The effect of interior density contrasts is 
to increase as deformation of internal boundaries supplements the 
perturbing potential* Compressibility decreases H^* 

Rather than solve for each Love number using a realistic earth 
model, I have chosen to use an approximate model calibrated using the 
known value of H 2 * The earth is assumed to be a sphere of 
incompressible fluid of density p with a point mass at the center 
sufficient to bring the average density of the sphere to p* The Love 
number for such a system is given by: 


H . ^ 

2£+l - 3p/p 


( 8 ) 


which reduces to (i + l/2)/(A - 1) for p ■ p, the value given by Munk 
and MacDonald for a uniform incompressible fluid* Using a value of 0*82 
for p/p gives the observed H 2 * This value, along with (7) was used to 
calculate the geoid anomalies used below. A more complicated two-layer 
model used to simulate bowing of Interior density constants, calibrated 
on H 2 , gives almost identical numerical results, although Its functional 
form is more complicated. 

As a result of the requirement that density contrasts be 
specified to second order, the thermal models used to describe the 
thermal evolution of the plates should be appropriate for plates which 
are spherical shells, rather than plane layers* The spherical analogue 
of the plate model consists of a spherical shell of thickness d which is 
initially at temperature T^. At times t > 0, the top of the shell (at r 
■ a) is held at T ■ 0, while the bottom (at r ■ a - d) is held at T^. 
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The temperature at depth z and time t Is: 

- ■'o T ♦ 7^ S «p(-n2.*«/d2)| (9) 


where < is thermal dlffusivlty (20). In the limit a * this becomes 
identical to the standard plate model with horizontal conduction of heat 
neglected (16). The steady>state temperature is less then, while the 
time-dependent temperature is greater than that for a plane layer plate. 

The corresponding expression for the boundary layer model, in 
which thn cooling penetrates arbitrarily deeply into the sphere is: 


T(z,t) ■ T^ 
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Only oceanic scafAoor less than 200 Myr old is of interest 
here, so all terms with tne exception of the positive, n > 0 term are 
negligible. Then 


T 


a - z 



( 11 ) 


The density contrasts leading to geoid anomalies in Eqs. (1) -(6) are 
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giver with respect to some arbitrary reference model* (The undulations 
of the geold are Insensitive to the details of the reference model)* I 
have taken an old, equilibrium thickness, oceanic plate as the standard 
for reference (Figure 1)* This model has 6.4 km of vater (p * Ig/cm'*) 
overlying a 120 km thick plate* The top 6*8 km of the plate Is crust 
with a density of 2*7 g/cm at 0°C; below the crust Is mantle of uniform 
composition* The teaq>erature varies through the plate according to the 
steady^atate part of Equation (9) from O^C at the top of the crust to T^ 
at the base of the plate, where the mantle density Is taken as 3*3 
g/cm^* Density Is related to temperature through a constant coefficient 
of thermal expansion* 

The density structures of the continental crust and 
lithosphere are not knotm, although variations can be bounded from 
observations of the geold (21). Beth vary with location* But since the 
primary purpose of this paper Is to conq>are models of the oceanic 
lithosphere, I have arbitrarily adopted a simple uniform continental 
structure which Is in isostatlc equilibrium with the reference oceanic 
model* 

Continental lithosphere Is assumed to be Identical to old 
oceanic lithosphere* Continental crust Is assumed to be 35 km thick and 
to extend 0*64 km above sea-level* For computational simplicity. Its 
density Is assumed constant at 2*7 g/cm^ In the upper 7.06 km, below 
which it has a decrease In density paralleling the variation with depth 
of the density ot oceanic crust caused by variation In temperature* 

For both the plate model and the boundary layer (sometimes 
called “cooling half -space") model, the crust is assumed to have 
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constant thickness. The depth of water Is constrained by the 
requlreaent that the lithosphere be In Isostatlc equilibrium with the 
equilibrium plate reference model. The density profiles for the 
three**type models are shown schematically In Figure 1. 

If equation (7) Is rewritten 
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then for the spherical analogue of the plate model, 
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For the boundary layer model. 
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In order to evaluate the integrala in (13) • (15) » It Is 
necessary to know the ago of all oceanic llthoaphere aa well at that 
part of the surface which is continental* Mauk (22) compiled the global 
eatiaatton of seafloor age used in this study on a 5 ** x 5^ grid* The 
expansions were carried through degree and order 20. 

Model geoids are obtained by summing the terms for each degree 
and order* Since geold anomalies are referenced relative to the center 
of mass of the earth, the degree zero and one terms ar ^3 not included in 
the sum* 

4* Isostatic Ceoid Models 

Although the primary purpose of this paper is to compare 
Isostatic geold anomalies for two models of the thermal structure of 
oceanic lithosphere, the computed geoids depend upon the model assumed 
for the subcontinental etructure* It is useful to Isolate the effect of 
this structure* Figure 2 shows the geold calculated using the 
continental structure in Figure 1 and neglecting the contribution from 
variation of oceanic lithosphere with age; in this model, all oceanic 
lithosphere is assumed to have reached thermal equilibrium. 

Continental regions in this model are characterized by geold 
highs of 5 to 10 meters while oceanic basins show lows of -5 to -10 m. 
The change in geold elevation of about 10 m at continental mergino 
matches that observed off the east coast of Australia (23), but is 
somewhat larger than that seen off the east coast of North Aaterlcan at 
40*S^N (8), indicating tliat the continental structure chosen is perhaps 
representative, but cannot be applied universally* The geold anomalies 
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in the North Atlantic are about 5 m more negative than those in the 
Pacific. This is a real effect of the three-dimensional nature of the 
distribution of continents as discussed later, not an artifact of the 
spherical harmonic expansion. 

The model geold which results when the "observed" variations 
in lithospheric age are used in the plate model (equation (13)], along 
with the standard continental model, are shown in Figure 3. Oceanic 
ridges and continents are characterized by geold highs of up to 6 m, 
while geold lows of about -7 m characterize the regions of Cretaceous 
seafloor in the Western Pacific, Atlantic, and Indian oceans. 

This Isostatic geold is similar to one computed by Chase and 
McNutt (24) using a spherical harmonic expansion of topography under the 
assumptions that continental topography ii compensated by the Airy 
mechanism at a depth of 35 km and that oceanic topography (excepting 
trenches) is compensated at 40 km. The main differences occur locally 
over high mountains, where their geold is about 10 m higher than mine, 
and trenches, where their geold is about 10 m lower. These differences 
are to be expected given the differences in the two reference models. 

The effect of variation in seafloor age alone, shown in Figure 
4, is obtained by subtracting the continental effect (Figure 2) from the 
combined oceanic age and continental model (Figure 3). The positive 
contributions from the midoceanic ridges become mc*e obvious. In 
addition, the contributions to the geoid highs from young lithosphere 
formed by back-arc spreading in the Western Pacific and off Sumatra 
become apparent. 

Figure 5 is the equivalent of Figure 3, this time calculated 
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using the spherical boundary layer model of the oceanic lithosphere* It 
Is remarkably similar to Figure 3, with the geoid anomaly over Jurassic 
oceanic lithosphere differing by only 2 m^ almost an order of magnitude 
less than the 15 m predicted on the basis of the one-dimensional models 
commonly used (17)* 

The marked similarity of the model geolds calculated using the 
full three-dimensional formulation applied to actual lithospheric age 
distributions is demonstrated directly by calculating the difference 
between the two models. Figure 6, which shows the difference between 
the boundary layer model (Fig* 5) and the plate model (Fig* 3), 
illustrates that the predictions of these two models is never much in 
excess of ± 2 m* Figure 6 employs a contour interval a factor of 5 
smaller than those in the previous models* 

5* Discussion 

The models usually used in discussing the Isostatlc geoid 
(7-9,17) are obtained by assuming 1) that the half-space equation (3) is 
applicable, and 2) that only local density contrasts (the k » 0 
wavenumber) is important* This is equivalent to assuming a 
one-dimensional density model* 

The spectra for some of the Isostatlc geoids shown above are 
presented in Figure 7* The solid line with solid dots gives the 
spectrum of the continental geoid (plotted in figure 2)* The spectra of 
the plate plus continent model (Figure 3) and boundary layer plus 
continent model (Figure 5) are similar to this one* The dashed line 
with open circles gives the net effect of the young lithosphere in the 
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plate model* (The geold for this model was shown In Figure A*) Both 
lines show a relative peak at 1 ~ 4->5* The spectrum for the difference 
between the plate and boundary layer models (Fig. 6) shown with 
crosses. Is about a factor of S smaller and doesn't reach Its relative 
peak until degree 10* Clearly the one-dimensional geold models commonly 
used are inadequate to describe geold models on a global scale, although 
they may be useful approximations near ridge-crests and continental 
margins. 

The plate and boundary layer models differ most for seafloor 
of age 80 My or greater. The largest contiguous area of the seafloor 
older than 80 My Is in the Western Pacific. This region has a diameter 
of roughly 35® of arc, which accounts for the peak in the spectrum of 
the difference between the plate and boundary layvir isostatlc geold 
models at 1 ~ 10. (The expansion Is carried through degree and order 
20. Thus, the small difference In predicted geolds Is a real effect, 
not an artifact of truncation of the spherical harmonic expansion). The 
very small difference between the theoretical geolds for the two thermal 
models arise because the long wavelength geold anomalies from ridges and 
continents common to both models extend Into the ocean basins and swamp 
the short wavelength differences between the models* 

The plate-driving force resulting from the thickening of the 
lithosphere with age Is proportional to the Integral of the product of 
the density contrast between the lithosphere and the asthenosphere and 
depth (25-27). As Parsons and Richter (17) pointed out, using a flat 
earth model, this driving force Is then proportional to the geold 
anomaly. Hager and O'Connell (27), using the plate model of the 
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lithosphere, quantified this drixing force on a spherical earth using 
Mauk's (22) distribution of lithospheric ages* They concluded that it 
!.s a significant driving force, distributed over plate area, which can 
help to explain the observed lack of dependence of oceanic plate 
velocity on area if plate motion is resisted by basal drag of a few 
bars* The calculations presented above show that their conclusions 
would not have changed significantly if they had used the boundary layer 
model rather than the plate model* 

Dahlen (18) advocates using as a definition of isostacy that 
1) shear stresses on vertical planes vanish and 2) devlatorlc stresses 
within the isostatically compensated layer be minimized* Requirement 2) 
leads to an increase in gravitational potential and geold height by a 
factor very close to (1 + 2)/t over that used here* I have chosen not 
to use this definition, preferring the definition used here which leads 
to a lower gravitational potential of the system* Minimizing devlatorlc 
stress would strengthen the conclusion of this paper, since it amplifies 
the low-degree harmonics more than the higher ones* It would also lead 
to a spectrum closer to that of the observed geoid* 

It is of Interest to compare the isostatlc geoids computed 
here with the observed geold heights* The top spectrum in Figure 7 is 
the spectrum of the observed geoid (28)* The observed spectrum is 
redder than the model spectra; at t « 3 the observed geold spectrum is 
a factor of 15 larger than the model ones, although at degree 20 the 
lithosphere contributes a large amount of power. 

Both thetmal models of the lithosphere have similar 
theoretical geoids when the actual distribution of seafloor ages is 
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used* Stripping the effects of seafloor aging using either model and 
the uniform continental reference model from the observed geoid should 
reveal the effects of dynamic processes In the mantle* Figure 8 shows 
the results of stripping the geoid in Figure 3 from the GEM8 geoid (28)* 
The GEMS geoid is referred to the hydrostatic figure, (f * 1/299*83), 
rather than the best-fitting ellipsoid, since any departure from 
hydrostatic equilibrium must be dynamically supported, and be of 
geophysical, if not geodetic, interest. 

The greatest amplitude anomaly in the stripped geoid remains 
over West Antarctica, with a value of -120 m* Central Siberia and 
Hudson Bay have anomalies in excess of -80 m* The low over Ceylon is 
barely -60 m, and appears to mark the edge of the anomaly stretching 
over the pole to the eastern U*S. The saddle in the geoid over Tibet is 
very likely the result of the doubling of crustal thickness there, which 
would lead to an additional local geoid elevation of about 20 m* 

The positions of positive geoid anomalies are not affected 
much, although there is a small decrease in amplitude over the Western 
Pacific (back-arc spreading removed) and Africa and S* America, where a 
small part of the highs are due to the effect of continental crust* 
Much (about 20 m) of the local geoid high over the Andes is the result 
of the great thickness there not included in the uniform crustal model* 
The main negative features of the residual geoid are two large 
anomalies; one is centered over W* Antarctica with an amplitude of 
almost -120 m, the other is larger in area, reaching from Ceylon, 
through Siberia to Quebec* There are also two large positive anomalies 
- one over New Guinea, the other centered just west of Africa* It has 
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been noted that these highs correlate well with subduct Ion and hotspot 
actlvltyt respectively (29,30)* The location of the Mt. Erebus hotspot 
near the center of the W* Antarctic low, the most pronounced geoldal 
feature, provides a noted and embarasslng exception to the correlation 
between geold highs and hotspots* 

The magnitude of the observed geold fluctuation Is an order of 
magnitude larger than those In the model geolds* Nonetheless, recent 
studies have shown that It Is possible to extract the lithospheric 
signal from the global "noise" In some areas (9,31). In the most 
comprehensive global study to date, Sandwell and Schubert (31) have 
shown that the slope of the geold height vs* age curve can be extracted 
In several regions* Near ridges these slopes are close to those 
predicted by one-dlmenslonal models using the plate or boundary layer 
thermal structures* Further from the ridges the slopes vary and even 
change signs* Although this type behavior Is not predicted by 
one-dlmenslonal models. Inspection of Figures 3 and 6 show 1) that this 
behavior Is to be expected given the three-dimensional nature of the 
problem, and 2) It cannot be used, given the "noise" In the geold 
signal, to discriminate between the plate and boundary layer models* 

Although plate and boundary layer models cannot be 
distinguished using global observations of the geold. It may still be 
possible to use geoldal variations across features such as fraction 
zones to place some constraint on the thermal structure of old 
lithosphere (32)* 
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Figure 1. 


Figure 2* 


Figure 3* 


Figure 4. 


Figure Ceptlons 


Density models used In calculating Isostatlc geold 
anomalies. A mature oceanic plate 120 km thick (right 
center and dashed lines) Is taken as a reference model* 
Continental structure (far right) Is assumed uniform; the 
subcontlnei'.tal lithosphere Is assumed to be equivalent to 
old ocean! r llthophere. Near the ridge (left center) the 
plate and boundary layer models are equivalent* For old 
lithosphere, cooling penetrates deeper Into the mantle for 
the boundary layer model (far left) than for the plate 
model* The Increased thickness of the lithosphere Is 
compensated by additional subsidence* 

Model Isostatlc geold resulting from the contrast In density 
structures between continents and oceans* The effects of 
oceanic age variation are not Included; all oceanic 
lithosphere Is assumed to be mature plate In thersial steady 
state* The contour Interval Is 5 m* 

Isostatlc geold resulting from the variation In thi^ oal 
structure with age of the plate model and from the contrast 
In density structure between continents and oceans* 

Isostatlc geold resulting from only the variation In thermal 
structure with age for the plat^ model* This figure 
represents the difference between the geolds In Figure 3 and 
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Figure 2* 

Figure 5* Isoetetlc geold reeulcing from the vavietlor. in thermal 

structure uxth age of the boundary layer model and from 
continents* 

Figure 6* Difference in isostatlc geoids for the boundary layer and 

plate models* obtained by subtracting the geold in Figure 3 
from that in Figure 5* The contour interval is only 1 m* 
If one-dimensional models were applicable, the maximum 
difference would be 15 m* 

Figure 7* Spectra for different geoids. The top line is the observed 

nonhydrostatic geold* The solid line with solid dots shows 
the effect of continental structure alone* The dashed line 
with open circles is the spectrum of the contribution of 
variations in oceanic age assumiPa the plate model* The 
bottom line represents the spectrum of the difference 
between the geoids of the bounder^ layer and plate models* 

Figure 8* Residual geold through degree and order 20 obtained by 

subtractng the isostatlc geold in Figure 3 from the GEMS 
geold referred to the hydrostatic flattening (f > 1/299*63)* 



















